;;########################################################################
;; mamdsmain.lsp
;; Multidimensional Scaling
;; Performs Metric MDS of a symmetricized dissimilarity data matrix or of a
;; Symmetricized Root-Mean-Squared Dissimilarity Data formed from several matrices
;; Copyright (c) 1991, 92 by Mary M. McFarlane
;; Written by Mary McFarlane, September, 1991
;;########################################################################

;;------------------------------------------------------------------------
;;function for constructing a multidimensional scaling plugin model object
;;------------------------------------------------------------------------

  (defun metric-averaged-mds 
    (&key 
     (data *current-data*)
     (dialog nil)
     (dimensions 3))
"Performs Metric Averaged Multidimensional Scaling.  With no arguments, performs a 3-dimensional scaling using all active data matrices in the current data. When there is more than one active data matrix, a single matrix is formed for analysis. It is formed by calculating, across matrices, the root-mean-squared value for every element in the data matrix. The result is made symmetric, if necesssary.
Keyword arguments are:
:DIMENSIONS followed by an integer to specify the dimensionality of the analysis (default - and minimum - is 3);
:DATA followed by the name of the data to be analyzed (default: *current-data*);
:DIALOG followed by t (to display parameters dialog box) or nil (default)."
    (if (not (eq current-object data)) (setcd data))
    (send cmds-model :new 
          "Metric Averaged MDS" data dialog dimensions))


(defproto cmds-model 
  '( dissim-object-name matrix-number matrices-numbers title data
         number-dimensions constant xmatrix old-matrix prev-matrix
         orig-matrix gamma transformed-data stress add-constants
         mds-2d-plot mds-3d-scatterplot mds-3d-spin-plot name-plot
         print-cmds-model stimulus-names subject-id vmenu) 
                () vista-analysis-plugin-object-proto)


(defmeth cmds-model :isnew 
  (analyze-menu-item-name data-object dialog dimensions) ;? MATRICES-NUMBERS
    (if dialog (setf dimensions (mds-options-dialog)))
    (when (or (not dimensions) (< dimensions 3))
          (send *toolbox* :reset-button button-number)
          (error-message "Minimum dimensionality is 3."))
  (setf dialog nil)
  (send self :matrix-number nil)
  (send self :matrices-numbers 
        ($position (send current-data :active-matrices '(all))
                   (send current-data :matrices)))
  (send self :add-constants t)
  (send self :number-dimensions dimensions)
  (send self :data-object data-object)
  (send self :mds-2d-plot-1 t)
  (send self :mds-3d-spin-plot-1 t)
  (send self :mds-3d-scatterplot-1 t)
  (send self :print-cmds-model-1 t)
  (send self :name-plot-1 t)
  (send self :fill-slots)
  (send (send *vista* :var-window-object) :clear)
  (send (send *vista* :obs-window-object) :clear)
  (send self :menu (send current-data :name))
  (send self :matrices (send current-data :active-matrices '(all)))
  (send self :dialog nil)
  (send self :icon-number (1+ (length (send *workmap* :icon-list))))
  (call-next-method analyze-menu-item-name data-object nil);dialog
  (send self :dissim-object-name (send self :name) )
  (send self :icon-title (send self :full-name))
  )

(defmeth cmds-model :analysis ())

(defmeth cmds-model :title (&optional (title nil set))
"Message args:  (&optional (title nil set))
It requires as its argument a string which it returns as the object's title."
  (if set (setf (slot-value 'title) title))
  (slot-value 'title))

(defmeth cmds-model :data (&optional (data nil set))
"Message args:  (&optional (data nil set))
It requires as its argument a list.  It returns a list of the data."
  (if set (setf (slot-value 'data) data))
  (slot-value 'data))

(defmeth cmds-model :stimulus-names (&optional (names nil set))
  (if set (setf (slot-value 'stimulus-names) names))
  (slot-value 'stimulus-names)) 

(defmeth cmds-model :stress (&optional (names nil set))
  (if set (setf (slot-value 'stress) names))
  (slot-value 'stress)) 

(defmeth cmds-model :add-constants (&optional (names nil set))
  (if set (setf (slot-value 'add-constants) names))
  (slot-value 'add-constants)) 

(defmeth cmds-model :number-dimensions (&optional (names nil set))
  (if set (setf (slot-value 'number-dimensions) names))
  (slot-value 'number-dimensions)) 

(defmeth cmds-model :mds-2d-plot-1 (&optional (names nil set))
  (if set (setf (slot-value 'mds-2d-plot) names))
  (slot-value 'mds-2d-plot)) 

(defmeth cmds-model :mds-3d-scatterplot-1 (&optional (names nil set))
  (if set (setf (slot-value 'mds-3d-scatterplot) names))
  (slot-value 'mds-3d-scatterplot)) 

(defmeth cmds-model :mds-3d-spin-plot-1 (&optional (names nil set))
  (if set (setf (slot-value 'mds-3d-spin-plot) names))
  (slot-value 'mds-3d-spin-plot)) 

(defmeth cmds-model :print-cmds-model-1 (&optional (names nil set))
  (if set (setf (slot-value 'print-cmds-model) names))
  (slot-value 'print-cmds-model)) 

(defmeth cmds-model :matrices-numbers (&optional (names nil set))
  (if set (setf (slot-value 'matrices-numbers) names))
  (slot-value 'matrices-numbers)) 

(defmeth cmds-model :matrix-number (&optional (names nil set))
  (if set (setf (slot-value 'matrix-number) names))
  (slot-value 'matrix-number)) 

(defmeth cmds-model :constant (&optional (names nil set))
  (if set (setf (slot-value 'constant) names))
  (slot-value 'constant))

(defmeth cmds-model :xmatrix (&optional (names nil set))
  (if set (setf (slot-value 'xmatrix) names))
  (slot-value 'xmatrix)) 

(defmeth cmds-model :old-matrix (&optional (names nil set))
  (if set (setf (slot-value 'old-matrix) names))
  (slot-value 'old-matrix)) 

(defmeth cmds-model :prev-matrix (&optional (names nil set))
  (if set (setf (slot-value 'prev-matrix) names))
  (slot-value 'prev-matrix)) 

(defmeth cmds-model :orig-matrix (&optional (names nil set))
  (if set (setf (slot-value 'orig-matrix) names))
  (slot-value 'orig-matrix)) 

(defmeth cmds-model :gamma (&optional (names nil set))
  (if set (setf (slot-value 'gamma) names))
  (slot-value 'gamma)) 

(defmeth cmds-model :dissim-object-name (&optional (names nil set))
  (if set (setf (slot-value 'dissim-object-name) names))
  (slot-value 'dissim-object-name)) 

(defmeth cmds-model :transformed-data (&optional (names nil set))
  (if set (setf (slot-value 'transformed-data) names))
  (slot-value 'transformed-data)) 

(defmeth cmds-model :subject-id (&optional (names nil set))
  (if set (setf (slot-value 'subject-id) names))
  (slot-value 'subject-id))

(defmeth cmds-model :name-plot-1 (&optional (names nil set))
  (if set (setf (slot-value 'name-plot) names))
  (slot-value 'name-plot))

(defmeth cmds-model :menu   (&optional (names nil set))
  (if set (setf (slot-value 'vmenu)  names))
  (slot-value 'vmenu))

;;______________________________________________________________________
;;given a (square) matrix, this function will calculate the euclidean
;;distances between each of the stimuli
;;______________________________________________________________________
(defun make-distances (x)
    (let* ((m (array-dimension x 0))
           (n (array-dimension x 1))
           (iden (identity-matrix m))
           (sum-prod (matrix (list m m) (repeat 0 (^ m 2)))))
    (dotimes (j n)
             (setf prod (outer-product (col x j) (col x j)))
             (setf prod-prod (* -2 prod))
             (setf di-mat (matrix (list m m) (repeat (diagonal prod) m)))
             (setf sum-prod (+ sum-prod 
                               (+ prod-prod (+ di-mat (transpose di-mat))))))
     (sqrt sum-prod)))
             
           
;;______________________________________________________________________
;;given a square matrix, this function will form the necessary v- and
;;v-plus matrices (Heiser and de Leeuw 1979, p. 9)
;;______________________________________________________________________

(defun make-vees (n)
  (let* ((iden (identity-matrix n))
         (one-matrix (- (+ 1 iden) iden))
    (v (* 2 (* n (- iden one-matrix))))
    (v-plus (* (/ 1 (* 2 n)) (- iden one-matrix))))
    (list v v-plus)))
    

;;______________________________________________________________________
;;given the square matrix of dissimilarities and a configuration matrix X, 
;;this function will make the
;;necessary B matrix found in the Heiser and de Leeuw (1979) article.
;;______________________________________________________________________

(defun make-bee (x disdata)
  (let* ((dist-matrix (make-distances x))
         (n (array-dimension x 0))
         (d-matrix (+ dist-matrix (identity-matrix n)))
         (b-matrix (* -2 (/ disdata d-matrix))))
    (dotimes (i n)
             (dotimes (j n)
                  (if (= i j)
                  (setf (aref b-matrix i j) (* -1 (sum (row b-matrix i)))))))
b-matrix))    
;;______________________________________________________________________
;;given a configuration x, this function will find the stress 
;;______________________________________________________________________


(defun calc-stress (x disdata)
  (let* ((dists (make-distances x))
         (num (/ (sum (square (- disdata dists))) (sum (square disdata))))
         (result (sqrt num)))
    result))

;;______________________________________________________________________
;;given a configuration x, this function will get the guttman-transform
;;______________________________________________________________________

(defun guttman-transform (x)
  (let ((n (array-dimension x 0)))
    (* (/ 1 (* 2 n)) (matmult (make-bee x (send current-model :transformed-data)) x))))


(defun square (x)
  (^ x 2))


  (defun add-constants-function (input-matrix)
    (let ((constant 0)
          (n (array-dimension input-matrix 0)))
        (dotimes (x (- n 2))
                 (dotimes (y (- n x 2))
                        (dotimes (z (- n x y 2))
                                 (let* ((i (aref input-matrix x (+ x y 1)))
                                        (j (aref input-matrix (+ x y 1) (+ x y z 2)))
                                        (k (aref input-matrix x (+ x y z 2))))
                                   (when
                                    (or 
                                     (< constant (- i j k))
                                     (< constant (- j k i))
                                     (< constant (- k i j)))
                                    (setf constant (max
                                                     (- i j k)
                                                     (- j k i)
                                                     (- k i j)))))))) constant))


(defun zero-diag (input-matrix)
  (let ((n (array-dimension input-matrix 0)))
    (dotimes (i n)
             (setf (aref input-matrix i i) 0))) 
  (/ input-matrix (standard-deviation input-matrix)))
  

(defun make-symmetric (x)
  (let* ((n (array-dimension x 0))
         (matrix-squared (^ x 2))
         (symmetric-matrix-squared (matrix (list n n) (iseq 1 (^ n 2)))))
    (dotimes (i n)
             (dotimes (j n)
                      (setf (aref symmetric-matrix-squared i j)
                            (/ (+ (aref matrix-squared i j)
                                  (aref matrix-squared j i)) 2))))
    (sqrt symmetric-matrix-squared)))


(defun find-bstar (x)
  (let* ((n (array-dimension x 0))
        (a2 (^ x 2))
        (means (/ (matmult a2 (rseq 1 1 n)) n))
        (matrix-means (outer-product means (rseq 1 1 n))))
        (* -0.5 (+ (- a2 matrix-means (transpose matrix-means))
                          (mean means)))))

(defun make-coords (bstar)
  (let* ((svd (sv-decomp bstar))
         (gamma (combine (select svd 1))))

  (let* ((sqrtgamma (sqrt gamma))
         (yoo (column-list (select svd 0)))
         (n (array-dimension bstar 0))
         (yoohoo (matrix (list n n) (iseq 1 (^ n 2)))))
    (dotimes (i n)
             (dotimes (j n)
                      (setf (aref yoohoo i j) (select (select yoo j) i))))
    (let* ((sqrtgammadiag (diagonal sqrtgamma))
           (xmatrix1 (matmult yoohoo sqrtgammadiag))
           (a (array-dimension xmatrix1 0))
           (b (array-dimension xmatrix1 1))
           (xmatrix (transpose (matrix (list (- b 1) a) 
                               (combine (select (column-list xmatrix1) 
                                                (iseq 0 (- b 2))))))))
(list xmatrix gamma)))))

(defmeth cmds-model :initialize ()
  (let* ((current-model-menu-length (length (send *model-menu* :items)))
         )
    (send *desktop* :no-menu-marks *model-menu*)
    (setf model-menu-item-name 
         (concatenate 'string  "MDS-" (send current-data :name)))
    (send *model-menu* :append-items
          (send menu-item-proto :new model-menu-item-name :mark t)
          )
    ))

(defmeth cmds-model :undo-iterations ()
  (let* ((n (array-dimension (send current-model :xmatrix) 0))
         (m (array-dimension (send current-model :xmatrix) 1))
         (prev (matrix (list n m) (copy-list (combine (send current-model :prev-matrix)))))
         (news (calc-stress prev (send current-model :transformed-data))))
    (send current-model :old-matrix 
               (matrix (list n m) (copy-list (combine (send current-model :xmatrix)))))
    (send current-model :xmatrix prev)
    (send current-model :new-sm-points)
    (send current-model :new-pp-points)
    (send current-model :new-sp-points)
    (send current-model :stress news)
    (send stressp :add-points
          (list (+ 1 (send stressp :num-points)))
          (list news) :color 'red)
    (send stressp :point-label (- (send stressp :num-points) 1)
          (format nil "~5,4f" news))
    (send stressp :add-lines-with-points :color 'red)
    (send stressp :point-state (- (send stressp :num-points) 2) 'normal)
    (send stressp :point-state (- (send stressp :num-points) 1) 'selected)
    (send stressp :redraw)))

(defmeth cmds-model :back-just-one ()
  (let* ((n (array-dimension (send current-model :xmatrix) 0))
         (m (array-dimension (send current-model :xmatrix) 1))
         (old (matrix (list n m) (copy-list (combine (send current-model :old-matrix)))))
         (news (calc-stress old (send current-model :transformed-data))))
    (send current-model :old-matrix (matrix (list n m) (copy-list (combine (send current-model :xmatrix)))))
    (send current-model :xmatrix old)
    (send current-model :new-sm-points)
    (send current-model :new-pp-points)
    (send current-model :new-sp-points)
    (send current-model :stress news)
    (send stressp :add-points
          (list (+ 1 (send stressp :num-points)))
          (list news) :color 'red)
    (send stressp :point-label (- (send stressp :num-points) 1)
          (format nil "~5,4f" news))
    (send stressp :add-lines-with-points :color 'red)
    (send stressp :point-state (- (send stressp :num-points) 2) 'normal)
    (send stressp :point-state (- (send stressp :num-points) 1) 'selected)
    (send stressp :redraw)))

(defmeth cmds-model :back-to-start ()
  (let* ((n (array-dimension (send current-model :xmatrix) 0))
         (m (array-dimension (send current-model :xmatrix) 1))
         (orig (matrix (list n m) (copy-list (combine (send current-model :orig-matrix)))))
         (news (calc-stress orig (send current-model :transformed-data))))
    (send current-model :old-matrix (matrix (list n m) (copy-list (combine (send current-model :xmatrix)))))
    (send current-model :xmatrix orig)
    (send current-model :new-sm-points)
    (send current-model :new-pp-points)
    (send current-model :new-sp-points)
    (send current-model :stress news)
    (send stressp :add-points
          (list (+ 1 (send stressp :num-points)))
          (list news) :color 'red)
    (send stressp :point-label (- (send stressp :num-points) 1)
          (format nil "~5,4f" news))
    (send stressp :add-lines-with-points :color 'red)
    (send stressp :point-state (- (send stressp :num-points) 2) 'normal)
    (send stressp :point-state (- (send stressp :num-points) 1) 'selected)
    (send stressp :redraw)))

(defmeth cmds-model :new-sm-points ()
  (let* ((n (send sm :num-points))
      (symbols (send sm :point-symbol (iseq n)))
      (states (send sm :point-state (iseq n)))
      (colors (send sm :point-color (iseq n))))     
    (send sm :clear :draw nil)
    (send sm :add-points  (col (send current-model :xmatrix) 
                               (iseq (send current-model :number-dimensions)))
          :point-labels (send current-model :stimulus-names)     
          :draw nil)
    (send sm :adjust-to-data :draw nil)
    (send sm :point-state (iseq n) states :draw nil)
    (send sm :point-symbol (iseq n) symbols :draw nil)
    (send sm :point-color (iseq n) colors :draw nil)
    (send sm :redraw)))

(defmeth cmds-model :new-sp-points ()
  (let* ((n (send sp :num-points))
      (symbols (send sp :point-symbol (iseq n)))
      (states (send sp :point-state (iseq n)))
      (colors (send sp :point-color (iseq n))))     
    (send sp :clear-points :draw nil)
    (send sp :add-points  (col (send current-model :xmatrix) 
                               (iseq (send current-model :number-dimensions)))
          :point-labels (send current-model :stimulus-names)     
          :draw nil)
   ; (send sp :adjust-to-data :draw nil)
    (send sp :point-state (iseq n) states :draw nil)
    (send sp :point-symbol (iseq n) symbols :draw nil)
    (send sp :point-color (iseq n) colors :draw nil)
    (send sp :redraw)))

(defmeth cmds-model :new-pp-points ()
  (let* ((n (send pp :num-points))
      (symbols (send pp :point-symbol (iseq n)))
      (states (send pp :point-state (iseq n)))
      (colors (send pp :point-color (iseq n))))
    (send pp :clear :draw nil)
    (send pp :add-points  (col (send current-model :xmatrix) 
                               (iseq (send current-model :number-dimensions)))
          :point-labels (send current-model :stimulus-names)     
          :draw nil)
    
    (send pp :point-state (iseq n) states :draw nil)
    (send pp :point-symbol (iseq n) symbols :draw nil)
    (send pp :point-color (iseq n) colors :draw nil)
    (send pp :redraw)
    ))

;  (LISTENERS)

(defmeth cmds-model :fill-slots ()
  (let* ((dobj (send self :data-object) )
         (temp (send self :transformer 
               (zero-diag (make-symmetric 
                           (if (eq nil (send self :matrices-numbers))
                               (send dobj :get-matrix
                                     (send self :matrix-number))
                               (/ (apply #'+ 
                                         (send dobj :get-matrices
                                     (send self :matrices-numbers)))
                                  (length (send self :matrices-numbers))))))
               (send self :add-constants)
               (send self :number-dimensions))))
    (send self :constant (select temp 2))
    (let* ((temp1 (make-coords (find-bstar (select temp 0))))
          (n (send self :number-dimensions))
           (m (array-dimension (first temp1) 0)))
      (send self :xmatrix 
            (transpose (matrix (list n m) (combine (select 
                     (column-list (first temp1)) (iseq n))))))
      (send self :prev-matrix (send self :xmatrix))
          ;  (transpose (matrix (list n m) (combine (select 
          ;           (column-list (first temp1)) (iseq n)))))
      (send self :old-matrix (send self :xmatrix))
      (send self :orig-matrix (send self :xmatrix))
           ; (transpose (matrix (list n m) (copy-list (combine (select 
            ;         (column-list (first temp1)) (iseq n))))))
      (send self :gamma (second temp1))
      (setf disdata (second temp)) 
    (send self :transformed-data (second temp)) 
    (send self :data (first temp)))
  (send self :stimulus-names 
        (send dobj :variables))
  (send self :subject-id (send dobj :matrices))
    (send self :stress (calc-stress (send self :xmatrix) 
                                    (send self :transformed-data)))))

(defmeth cmds-model :transformer (use-matrix add-constants number-dimensions)
     (setf constant (add-constants-function use-matrix))
     (setf trans-matrix (zero-diag (+ use-matrix constant)))
     (list use-matrix trans-matrix constant))

; (LISTENERS)

(defmeth cmds-model :options ()
  (when (send self :dialog) (mds-options-dialog)))

(defun mds-options-dialog ()  
  (let* ((d-label (send text-item-proto :new "Dimensions (minimum=3)"))
         (dimensionality-text
          (send edit-text-item-proto :new "3" :text-length 1))
         (cancel (send modal-button-proto :new "Cancel"))
         (ok (send modal-button-proto :new "OK" 
                   :action #'(lambda ()
                       (eval (read (make-string-input-stream
                                          (send dimensionality-text :text))
                                         nil)))
                   ))
         (mds-options-dialog 
          (send modal-dialog-proto :new
                (list (list dimensionality-text d-label)
                      (list ok cancel))
                :default-button ok))
         )
    (send mds-options-dialog :modal-dialog)
    ))

(defmeth cmds-model :save-model-template (data-object)
"Args: (data-object)
DATA-OBJECT is the object-identification information of a data object. 
The method contains a template for saving the model-object."
  `(multidimensional-scaling    
    :title      ,(send self :title)
    :name       ,(send self :name) 
    :dialog       nil
    :dimensions ,(send self :number-dimensions)
    :data (data ,(strcat "Analyzed " (send data-object :name))
                :title      ,(strcat "Analyzed " 
                                     (send data-object :title))
                :variables ',(send self :stimulus-names)
                :labels    ',(send self :stimulus-names)
                :matrices  ' ("Analyzed Dissimilarity Data")
                :data      ',(combine (send self :data))))
  )

(defmeth cmds-model :report (&key (dialog nil))
  (if (not (eq current-object self)) (setcm self))
  (let* ((w nil)
         (stress (send self :stress))
         (constant (send self :constant))
         (eigenvalues (send self :gamma))
         (proportions (/ eigenvalues (sum eigenvalues)))   
         (stimulus-names (send self :stimulus-names))
         (nstim (length stimulus-names))
         (ndim (send current-model :number-dimensions))
         (fitmat (transpose (matrix (list 3 nstim) 
                                    (combine eigenvalues proportions
                                             (cumsum proportions)))))
         )
    (setf w (report-header (send self :title)))
    (display-string (format nil "Multidimensional Scaling Analysis~2%Model: ~a" (send self :name)) w)
    (display-string (format nil "~2%Variable (Stimulus) Names: ~a~%~%" stimulus-names) w)
    (display-string (format nil "Analyzed Dissimilarity Data (with Additive Constant).~%Derived from") w)
    (if (= 1 (length (send current-model :matrices-numbers)))
        (display-string (format nil " Matrix ~a~%" (send self :matrices)) w)
        (display-string (format nil " the Average of Matrices ~a~%" (send self :matrices)) w)
        )
    (print-matrix-to-window (fuzz (send self :data)) w :labels stimulus-names)
    (display-string (format nil "~%Additive Constant: ~g~%" constant) w) 
    (display-string (format nil "~%Current Stress: ~g~%" stress) w) 
    (display-string (format nil "~%Eigenvalue (amount of data variance fit), Proportion & Cumulative Proportion of variance fit for each dimension~%") w)
    (print-matrix-to-window fitmat w)
    (display-string (format nil "~%Initial Stimulus Coordinates:~%") w)
    (print-matrix-to-window (fuzz (select (send self :orig-matrix) (iseq nstim) (iseq ndim)) 3) w :labels stimulus-names)
    (display-string (format nil "~%Current Stimulus Coordinates:~%") w)
    (print-matrix-to-window (fuzz (select (send self :xmatrix) (iseq nstim) (iseq ndim)) 3) w :labels stimulus-names)
    (send w :fit-window-to-text)
    w))

(load (strcat *ma-mds-plugin-path* "mamdsmob2.lsp"))